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We have numerically studied the statics and dynamics of a model three-dimensional vortex lattice 

at low magnetic fields. For the statics we use a frustrated 3D XY model on a stacked triangular 

lattice. We model the dynamics as a coupled network of overdamped resistively-shunted Josephson 

junctions with Langevin noise. At low fields, there is a weakly first-order phase transition, at which 

'j-y-T ' the vortex lattice melts into a line liquid. Phase coherence parallel to the field persists until a sharp 

Qs I crossover, conceivably a phase transition, near Ti > Tm which develops at the same temperature 

Q^ ■ as an infinite vortex tangle. The calculated flux flow resistivity in various geometries near T = Ti 

^— ( ' closely resembles experiment. The local density of field induced vortices increases sharply near 

J, , Te, corresponding to the experimentally observed magnetization jump. We discuss the nature of a 

C^ ■ possible transition or crossover at TeCB) which is distinct from flux lattice melting. 
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^^ ' Ever since their discovery, the behavior of high-Tc materials in a magnetic field has seemed mysterious |ll| . Unlike 

^\J ' the conventional low-Tc type-II materials, high-Tc superconductors (HTSC's) show a broad region in the magnetic- 
C\| ■ field/temperature (H-T) plane where the Abrikosov lattice has apparently melted into a liquid state M . 
^~~l ■ Considerable recent evidence now suggests that flux lattice (FL) melting is a first-order phase transition. On 

^^ the experimental side, a local magnetization jump has been measured by network of Hall microprobes H on 
Bi2SrCa2Cu208 and has been associated with the melting transition. The transition thus observed seems to lie quite 
near the melting curve as determined from low angle neutron diffraction [Q and /iSR experiments 1^. More recently, 
Schilling et al MM have directly observed the latent heat of the transition in YBa2Cu307„5 along a line T„i{H) in the 
H ^ T phase diagram which agrees well with mechanical and transport measurements [p|-|lO[ . Numerical evidence for 
i-rt ' a first order melting transition has been obtained from simulations based on a frustrated XY model [HI 12 1, and from 
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^ ' a lowest Landau level model which is expected to be most accurate at high magnetic field |13|. First order melting 
O ' has also been observed numerically in a system of unbreakable flux lines described by a Lawrence-Doniach model mm . 
^ ' All these simulations are based on a large density of flux lines {^ 0(1 — 10)tesla). 

^ , An anomalous feature of the local Hall probe measurements is that the apparent first order transition line seems 

to terminate at a critical point above which the latent heat vanishes [p|. Since on symmetry grounds a first order 



^ . , "melting" line cannot terminate in a critical point [hSl , this critical point may suggest that the first order melting line 
j^ ■ is instead intersecting another phase transition line related to the disorder. A related issue is the entropy released per 

vortex per layer across the transition line. This entropy increases very rapidly as the field decreases. Such behavior 

is difficult to accounted for within a model based only on the field induced vortices. 

In the presence of disorder, the lattice becomes unstable at high fields against proliferation of quenched-in topological 

defects [|6 17 , possibly through a first order phase transition across a horizontal (constant H) line in the H-T plane. 



This line then may meet the temperature-driven melting line, causing it to terminate. Somewhere along this line, 
the melting transition may be converted into the universality class of the continuous vortex glass transition p8[ , 
characterized by divergent correlation lengths and times. 

Another unresolved issue regarding the phase diagram is the possibility of reentrant melting at low fields. Reentrant 
flux lattice melting is expected because of screening of the widely-separated vortex lines at low flelds ||,|l^. It has 
been recently reported in single-crystal NbSe2 sample [p^] , based on tracking of the so-called "peak effect" pl|,p2| . 
Such reentrance behavior has been observed only at the limited field range by Ling et al [ p3[ . On the other hand, the 
melting curve tracked by the micro-Hall probe ^J seems to monotonically approach the zero-field superconducting 
transition at Tc{H = 0) even for fields as low as a few Gauss. 

FL melting can also be probed by transport measurement. But since such measurements are non-equilibrium, they 
offer only an indirect means of studying equilibrium FL melting. In real materials with disorder, the interpretation of 
transport measurement is further complicated by the many competing energy scales. In single-crystal YBa2Cu307_5 , 



the in-plane resistivity exhibits a discontinuous jump and hysteresis which have been identified with a first-order 
mehing transition |8|. Nonetheless, the peak effect in the critical current occurs at slightly lower temperatures than 
the resistivity jump, leading some workers to postulate that there is a "premelting" phenomenon |2^] in this material, 
in addition to melting. In Bi2SrCa2Cu208 , simultaneous transport and local magnetization measurements ^4| show 
that the jump in local magnetization M coincides (at high fields) with a jump in the resistivity pab from zero to a 
finite value, or (in low fields) the continuous development of a finite pab- In addition, at high fields, the jumps in p 
and M are accompanied by hysteresis. Together, these phenomena strongly suggest first order flux lattice melting at 
high fields. At low fields, the experiments are more ambiguous. 

FL melting has been widely studied numerically. The possibility of two stage melting was first suggested by Li and 
Teitel [£5[g6| for a model with infinite penetration depth A, and later for a system with finite A |2^,^. The calculations 
of Li and Teitel are based on the so-called frustrated XY model with fairly low flux per plaquette of / = 1/25 (in 
units of the flux quantum $0 = hc/2e) on a simple cubic lattice. They find that the three-dimensional flux line 
lattice (3D FLL) melts first into a "line liquid" characterized by disentangled flux lines, which become entangled at a 
second, higher-T phase transition. Current-voltage (IV) measurements in the so-called "flux transformer" geometry 
|P9|-p2| provide some support for this picture. Specifically, they suggest that FL melting is signaled by the onset of 
finite in-plane resistance, while in an applied current, phase coherence is lost in the c direction only at a distinctly 
higher temperature. On the other hand, simulations of dense (f=l/6) fiux lines on a stacked triangular grid favor a 
single transition ||l|,|l^. Dynamical calculations |g^ on a triangular lattice at / = 1/6 suggest that if there are two 
separate transitions, they arise from pins, either intrinsic to the discrete cubic grid, or put in by hand. Yet more recent 
studies based on a London vortex loop model on a simple cubic lattice show that superconducting order disappears 
apparently in two steps, the sequence of which depends on the lattice anisotropy |p4|, although it is argued to be a 
finite size effect by the same authors p5| ]. 

In this paper, we attempt to resolve some of these issues by considering the frustrated XY model over a wide 
range oj flux densities, using both static and dynamic simulations but with no quenched disorder. By examining this 
model on a stacked triangular lattice, we minimize the unphysical periodic pinning due to the lattice. By working at 
relatively low densities, we focus on the regime, now being probed experimentally, where the XY phase fluctuations 
(vortex loops) are as important as those oi field induced vortex lines [ p6[ . Our main conclusion is that there are, 
in fact signatures of two separate transitions at low fields, which are not artifacts of pinning by the discrete grid. 
The transition at lower temperature is unambiguously associated with vortex lattice melting. The second transition 
may be a sharp crossover rather than a true phase transition. Nevertheless, it is responsible for several experimental 
features (such as sharp increases in local magnetization and in resistance) which are often identified as evidence for a 
first order melting transition. 

The remainder of this paper is organized as follows. In Section O, we describe our model and its numerical solution. 
The following sections present our numerical results, which are followed by a discussion and then summarized in a 
concluding section. 

II. MODEL 
A. Hamiltonian and Thermodynamics 

We study the standard frustrated XY model described by the Hamiltonian 

W = -J^cos(0, -0,-yl„), (1) 

where Aij — ^ J^ A ■ dl, A is the vector potential associated with a uniform magnetic field B = Bz applied parallel 
to z, $0 = hc/2e is the flux quantum, 6i is the phase of the order parameter on site i, and the sum runs over nearest 
neighbor pairs. We use a stacked triangular grid with B||z, the direction perpendicular to the triangular network, 
with periodic boundary conditions (PBC) in all directions except where stated otherwise. 

To allow a wider range of frustrations compatible with the boundary conditions, we use a variant of the Landau 
gauge Ig^. Note that there are four bonds per grain: three in the xy-plane and one along z. We label these by 
their unit vectors a = x, 2/1,2/2, z where 7/1 = (l/2)i; -I- ('\/3/2)y and 2/2 = —{V3/2)x + {l/2)y. The phase factors Aij 
connecting a grain located at (x,y,z) to its four nearest neighbors are given by along x or z, 2TTf{2x+ 1/2) along 2/1, 
and 27r/(2a; — 1/2) along 2/2- There are exceptions to this form for grains lying on the boundaries: All grains lying 
on the X — hx boundary plane have Aij = — 27r • 2Lixy for bonds in the x direction. Bonds at a; = La; boundary such 
that mod(j, 2) = 1 have Aij = 2TTf[2x + 1/2 — 2Lx{y + 1)] in the 2/1 direction. For bonds on the x = boundary with 



mod(j, 2) = 0, Aij = 27r/[2x — 1/2 + 2Lx{y + 1)] for bonds in the 2/2 direction. In contrast to the usual Landau gauge 
(which is compatible with frustrations only in integer multiples of l/{2Nx)), this generalized gauge is compatible with 
any / which is an integer multiple of l/{2NxNy) under periodic boundary conditions. 

We have considered networks of sizes Af = N^ x Ny x N^. For / = 1/24, we have studied Nz = 12,24,48 and 
N^ = Ny = 24, and for other values of / (1/2592, 1/1648, 1/81, and 1/6) we have considered N.^ ^ Ny ^ N^ = 18. 
In two dimensions, the vortices lie on the vertices of a honeycomb grid of unit length of (1/V3)as which is dual 
to the triangular grid of unit length a^. Assuming that vortices form perfect triangular lattice on this grid, and 
equating the area per vortex to (•\/3/4)a^//, we obtain the following necessary condition for a triangular vortex 
lattice to form without geometric frustration of the FL:2// — (tii/3 + ri^/'i) with integers rii, n2. The values of 1// 
satisfying this condition are then 2, 6, 8, 14, 18, 24, 32, 38, 42 ... , 648, . . .. For / = 1/162, two distinct pairs [(ni, 712) 
= (0,36), and (27,18)] satisfy the condition. / — 1/648 has the lowest possible nominal vortex density of one per 
18 X 18 system compatible with our chosen gauge, and allows either of the pairs (71.1,712) = (0,72) and (54,36). 
/ — l/2592(= 1/4 X 2x18x18 )' represents less than a single vortex line, and the system is gauge-frustrated. 

In practice, for / < 1/81, there are too few field-induced flux lines to study FL melting. Nonetheless, the dilute 
regime is still of interest, since in these cases, the flux lines behave independently and the thermodynamics is dominated 
by the zero-field phase degrees of freedom p^. Except for f = 1/2592, we study only gauge-unfrustrated values allowing 
only an integer number of vortices in the simulation box. 

We calculate the thermodynamics using a standard Monte Carlo (MC) algorithm, with up to 10^ MC steps at each 
temperature T. To ensure equilibration in the ground state for all values of /, we performed simulated annealing runs 
for the two-dimensional (2D) version at each / with the same lateral dimensions. We then form the ground state of 
the 3D system by stacking the 2D ground states thus found uniformly along the z— direction. This enables us to find 
the ground state configuration of a perfect triangular lattice for low values of 1/24 < f < 1/18 within a reasonable 
time. Starting from these 3D ground states, we warm up the system in steps of AT/ J ~ 0.05 or 0.1, allowing at least 
4 — 5 X 10^ Monte Carlo sweeps for each T. The final configuration for each T is then saved to be used as a starting 
configuration in some of the longer calculations as well in the dynamic simulations. 

From these calculations, we extract a range of thermodynamic quantities. One of these is the specific heat Cy = 
{{H^) — {H)'^)/{kBT) at temperature T. We also calculate the local vorticity vector field na{p) deflned for each 
Cartesian direction a and each point p of the stacked honeycomb dual lattice. At each instant during the simulation, 
na{p) is determined from 
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Aij,2^]=27r[n„(p)-fp]. (2) 



Here the summation runs along the bonds belonging to the plaquette labeled p, a (a triangle in the xy plane, a 
square in planes parallel to the z axis); and fp = ^ Aij / {2t:). From tt-q, one can also compute the structure factor, 
S'a/3(k) = (7ic((k)r7^(— k)), where nQ,(k) is the Fourier transform of the local vorticity vector nct(r). We also calculate 
the principal components ^xx and ^zz of the helicity modulus tensor ||38| , in the directions perpendicular and parallel 
to the applied field. To within a constant factor, 7 represents the phase rigidity or the superfluid density tensor of 
the system; its derivation in terms of equilibrium thermodynamic averages has been given elsewhere [p9| . 

B. Dynamics 

To treat the dynamics, we model each link between grains as an overdamped resistively shunted Josephson junction 
(RSJ) with critical current I^ = 2eJ/h, shunt resistance i?, and Langevin white noise to simulate temperature 
effects. The effective IV characteristics are then obtained by numerically integrating the coupled RSJ equations, as 
described elsewhere pO[, using a time constant typically of O.lio and obtaining voltages by averaging over an interval 
of ~ 600io — 2000io- Since direct solution of these equations would involve inverting and storing an A/" x A/" matrix, 
where M is the total number of grains~ 0(5000), we instead solve them iteratively |4^, incurring a speed penalty of 
a factor of In A/". We verified that our solutions converge by comparing them with those from direct inversion for time 
steps of 0.01, 0.04 and 0.1 on an 8 x 8 x 8 system. 

The most obvious approach to the dynamics of this model would be to use free boundary conditions, injecting 
current into one face of the lattice and extracting it from the opposite, with periodic transverse boundary conditions 
J42[ . But this has the following disadvantage. Once the flux lattice is depinned from its underlying periodic pinning 
potential, it will drift along as a whole under the influence of the Lorentz force provided by the driving current. 
Since this occurs equally in the solid and the liquid state, such a geometry may not distinguish clearly between flux 
lattice and flux liquid (in the absence of spatially inhomogeneous pinning centers). This problem may be even more 



conspicuous in our stacked triangular geometry, since the critical current Idp for depinning a single vortex pancake 
from the underlying triangular grid at zero temperature (T = 0) is smaller {Idp ~ 0.037/c) than in a square grid 



(7rfp«0.1/,) |43|. 

We therefore adopt a different geometry for injecting and extracting current, as shown in Fig. ^. Fig. ^ (a) 
corresponds to injecting current I / Ic into each grain in the yz plane at x = 0, and extracting it from each grain at 
X — Nx/2 (with periodic boundary conditions in all three directions). In this geometry, the Lorentz forces acting on 
the vortices in the two halves of the volume are oppositely directed, as indicated by the arrows. Thus, in this geometry, 
we are effectively probing the shear modulus /i of the vortex lattice, on a length scale La;/2. Similar geometries have 
been previously discussed in the context of possible experiments Q,Q. In Fig. |l| (b), we show a geometry which is 
designed to probe the c-axis resistivity, pc- In this case, we inject a current / into each grain on the xy plane at z = 
and extract it from each grain at z = N^/l. There is on average no Lorentz force on the vortex lines. 

III. ZERO-FIELD XY MODEL: F = 
A. Thermodynamics 



The 3D unfrustrated XY model on a cubic grid has been extensively studied |46|. Near the phase transition, the 
specific heat Cy - |T - TxfT" with a ~ and the correlation length C ~ |T - Txy\~'' with v ~ 0.66 - 0.67. For a 
cubic lattice, Txy ~ 2.203 J. In a stacked triangular grid, where each grain has more nearest neighbors, the transition 
is shifted to a higher temperature. Numerically, we find that Txy '^ 3.04J. 

The XY phase transition is best characterized by the helicity modulus tensor j^fs, which measures the phase rigidity 
of the system |3^]. In stacked triangular lattice, this tensor is diagonal with elements 
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( (E sin(9y)?T,y ■ na]'^) ~ ('^ sm{Qij)fiij ■ ha) Y (3) 

Here A* is a fictitious uniform vector potential in the a direction, V is the volume, Qij = Oi — 9j — Aij is the 
gauge-invariant phase difference, n^j and ha are unit vectors along the ij*"^ bond and in the a direction. 

It is useful to distinguish two contributions to @ij'- one due to spin waves, and one due to vortices p7|. The 
former is dominant when the sinG^ « 0^^, while the other is nonzero when the vorticities na{p) ^ 0. Thus we write 
laa = laa^ +lXaj whcrc the two terms on the right hand side are respectively the spin wave and vortex contributions 
to 'Yaa- The spin- wave degrees should predominate at low temperatures, while the vortex degrees of freedom are the 
dominant excitations near the phase transition |4q-p2]. The spin- wave contribution 7'^^ can be estimated within a 
self-consistent harmonic approximation with the result |53| 



T^"* ~ Jexp 
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where Z? = 4 for a stacked triangular lattice. 

To characterize the vortex contribution, we introduce the net global vorticity vector by p4 

■Ma — ha ■ daOy 



e^da- / e^da (5) 

2+ JT,- 

where S+ and S^ are the two bounding planes normal to ha, with normal vectors parallel or anti-parallel to ha- 
We assume that the singular portion of the phase variable 9y has been selected out. A^^ is sensitive to existence of 
unbound vortex lines perpendicular to ha- This is illustrated in the left panel of Fig. g for the case of a single infinite 
vortex line piercing the sample normal to the a direction. In this case, the phase integrals on the planes S+ and S~ 
give nearly equal but opposite values. Thus Ma has large fluctuations, leading to a reduction in the value of ^aa (see 



below). Closed vortex loops, such as shown in the right panel, give a zero contribution to Ada- In general, for f — 0, 
the thermal average (Ma) ^ 0. For an applied field || z, (^ Alc\ ^ for a = a; or y. It can be shown that ^aa and 
Ma are related by 

1 J2 r/. ,.\ /. , \2i 



}aa laa 
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;{{Ml)-{Ma)}. (6) 



Thus, the vortices make a negative contribution to jaa arising from fluctuations in M. They may be said to 
predominate over the spin waves when their fractional contribution to the helicity modulus is of order unity, that is 



J 



{{Ml)-{May}/Af-^0{l) (7) 



where J\f is total number of grains. 

In Fig. 0, we show the calculated 7^^, as well as the value j^z^ determined from the self-consistent harmonic 
approximation (SCHA), eq. (j^. The other principal components of 7 behave similarly. 7^^ and 'y^^ begin to differ 
for temperatures as low as T ~ O.STxy, where vortex loops start to be excited. The SCHA predicts a discontinuous 
jump in "fsw from a value of about 0.37 at Txy to zero. This jump is an artifact of the approximation, which neglects 
the periodicity of the Hamiltonian in the angle variables and the vortex fluctuations. The inset shows a finite size 
scaling analysis to locate Txy- The helicity modulus 7 ~ \T — Txy\^ with v = {d—2)v. Therefore, the scaled quantity 
7L for an L X L X L system should cross a single point at Txy- Based on this criterion, our numerical results give 
Txy = 3.04J ± 0.02. We also observe that 7(r) approaches zero with v - 2/3 for 0.02 < \Txy - T\/Txy < 0.1 and 
deviates from this outside the range. 

Eq. (0) is equivalent to that derived in Fourier space by Chen and Teitel pq| for A —> 00: 



VT q 



if we take q —^ n/Jj^ ■ 

An alternative form based on the vortex loop scaling picture has been obtained by Williams pq | in the limit A — > 00 
where Ta^ is given in terms of vortex loop diameter distribution. To check the importance of large loops, we show 
in Fig. H, our calculated size distribution of connected vortex loops near Txy- Two vortex segments are considered 
connected if they cross within a single unit cell. Such crossings become very extensive at Txy, suggesting that the 
energy barrier for vortex line intersections vanishes near Txy For T < Txy, the maximum vortex loop size is finite, 
while for T > Txy, there start to occur vortex lines spanning the entire simulation cell. Thus Txy somewhat resembles 
a bond-percolation transition, although Txy does not correspond to the percolation threshold, but to a point where 
connected cluster first forms a D — 1 dimensional manifold if the system is in dimension D. Because such infinite 
clusters occur, it appears that the behavior seen in Fig. H is not a consequence of the finite simulation cell, but persists 
in the thermodynamic limit. A similar picture, but with no long-range interactions among the vortex segments (the 



so-called polymer limit), has been discussed by Akao |5q] and by Kultanov et al |57 



B. Dissipation near Txy 



We have also calculated the dissipation near TxY^f = 0), using the periodic current injection geometry of Fig. H, 
In this case, since there are no field-induced vortex lines, the calculated dissipation can be unambiguously related 
to the resistivity. Our results ar e show n in Fig. for several values of the bias current density. [More details of the 
method are described in section VI D| .] To calculate the differential resistance, dV/dl, we carried out two separate 



runs at I/Ic = 0.083 and 0.043, to obtain dV/dl = y(0.083) - y(0.043)/0.04. Fig. | shows this result as well as 
R = V {0.083) /I in the inset. These bias currents are low enough to show sharp features at TxyH = 0) while not 
significantly disrupting that transition. At higher current densities (not shown), there are numerous current-induced 
vortex loops. These increasingly round out the sharp jump at Txy if — 0) shown in the Figure, which eventually 
washes away entirely. 

Fig. H shows the average number of vortex segments per plaquette as calculated both by Monte Carlo simulations 
(with no driving current) and by coupled RSJ dynamics (with a finite bias current). Evidently, just at the XY 
transition, the system becomes filled with thermally generated vortex segments, one per grain, or elemental cell. Below 
TxyH = 0), while there may exist vortex loops of arbitrary size, the number of these, V{1), falls off exponentially 
for large clusters (cf. Fig. H). By contrast, for T > Txy, "^(0 diminishes algebraically with I. This subtle change in 



Vjl) implies that the average size of the connected vortex tangle remains finite for T < Txy, but diverges above it 
PSj . We find that there are numerous finite vortex loops at temperatures as low as 0.5Txy- Moreover, at any given 
temperature, we find that a large bias current further enhances both their number and their size. 

A simple argument suggests that the dissipation below TxyH = 0) may be exponentially activated. The ef- 
fective energy for a vortex loop of radius r oriented normal to a uniform driving current density j is U{r) ^ 
2'iTr ■ Min[ln(r), ln(A)] — c • Trr^j, where A is the penetration depth. Thus, there is a barrier to loop expansion 
with a critical radius re ^ In A/j and height Umax ^ (In A)^/j. For sufficiently small j, only those vortex loops of size 
r > re will expand and contribute to dissipation. From our simulations, for T < Txy, ^(0 decays exponentially. 
Therefore, the small dissipation involving the expansion of thermally nucleated vortex rings should have an activated 
temperature-dependence for T < Txy, resulting in highly nonlinear IV characteristics. 

IV. DENSE LIMIT: F = 1/6 

This is the largest value we studied which allows a triangular vortex lattice commensurate with the underlying 
triangular grid. This value yields a strong first order transition pi|jl^ ], with an entropy of melting AS* of about 0.3 
ks per vortex pancake. 

In the present MC simulation, the lattice was gradually warmed up from a perfect triangular lattice in temperature 
steps of dT — O.IJ, with 50,000 MC sweeps for each T. The insets show the in-plane density-density correlation 
< nz{r±, z)nz{0,0) > for T/J — 1.175 and 1.3, slightly below and slightly above the melting temperature. At the 
melting temperature T™, our MC histogram for the internal energy distribution agrees with that of ||ll[. We have 
also calculated both ^zz{T) and the Bragg intensity Szz(Gi) at the smallest reciprocal lattice vector Gi for the 
triangular vortex lattice. The results are shown in Fig. 0. To within 5T / J ^ 0.025, both quantities vanish close 
to T = Tm(l/6) ^ 1.175J, the melting temperature as determined from the double peaks in the energy histogram 
im . The apparent melting T k, 1.2 J, slightly higher than inferred from the energy histogram, seems to be due to a 
superheating effect. 

The occurrence of only a single phase transition at f — 1/6 [ pTpq | is not surprising: at this field, there is one vortex 
pancake for every three grains, and hence, only three phase degrees of freedom per pancake. Thus, most potential 
vortex excitations are already exhausted by the lateral fluctuations of field-induced vortex lines in the liquid phase. 
This can be seen in Fig. @, where we show two typical vortex configurations, one slightly below and the other above 
the melting transition. Clearly, the transverse line fluctuations quickly dominate the thermodynamics above Tm- 
The strong first order transition at / = 1/6 can be then understood by the close connection between the lateral line 
fluctuations and incipient vortex loops. We conclude, with an accuracy of dT/ J = 0.025, that in the dense limit of 
/ — 1/6, superconducting coherence is destroyed in all directions as soon as the lattice melts. 

V. DILUTE LIMIT 

By dilute limits we mean the regime where the number density of thermally excited vortex line segments Uc equals 
or exceed the density of field-induced vortex line segments. To make this more quantitative, we first consider a perfect 
line lattice at T = with a given /. The number of unit vortex line segments per grain will be 1/(3/), or 1/(15/) per 
plaquette (since there are five plaquettes per grain). Thus, for / = 1/6, there are 0.4 vortex segments per plaquette. 
In Fig. 0, we show the density Uc of thermally excited vortex segments at / = 0, including all three directions. Note 
that at the XY transition, Uc = 0.15. Thus, / = 1/6 is clearly in the dense regime, while / < 1/18 is roughly in the 
dilute regime. 

Fig. P shows the specific heat Cy per grain, as calculated from energy fluctuations for several values of f [1/162 (4 
flux lines), 1/81, 1/24, 1/18 and 1/6 (108 lines)] in both the dilute and dense regime. At low T, all the Cy's approach 
ks/i per grain, as expected from the Dulong-Petit law. The overall behavior of the peaks in Cy up to / = 1/24 is 
remarkably similar to that seen in YBa2Cu307_5 [q8[. For / = 0, it is known that Cy has a weak divergence (cxjtj__" 
with a « 0.0 iQ), as expected for the 3D XY model. At flnite /, this peak is rounded as seen experimentally p8| , p9| . 
As discussed below, these broad peaks in Cy generally occur well above the melting transition in this field range. 
Note that our results for the dense (/ = 1/6) case differ qualitatively from all those at lower /. The sharp peak 
for / = 1/6 is actually a delta- function singularity, consistent with the finite heat of fusion of a first-order transition 
known to occur at this density pJ^ . 



Fig. 10 shows that the height of the peak for / < 1/18 roughly follows a logarithmic dependence on the magnetic 
length Lb defined as Lb ~ l/VTi tbe average vortex spacing. Furthermore, as we show in Fig. |l^, the position of the 
peak at a finite / shifts from TxY{f = 0) by an amount STc{Lb) which closely follows the law ~ Lg'' . As we will 



show in more detail for / = 1/24 below the phase rigidity along the applied field, as measured by 7^2, vanishes for all 
/ near the broad maximum in specific heat. At the temperatures [which we denote Ti{f)] where jzzif) vanishes, we 
observe that the average number of thermally generated vortex segments per grain normal to the z direction closely 
follows the law n'^ (Tg) ^ 0.1 • L^^*°'^. All this behavior is discussed in more detail below. 

VI. F=l/24 
A. Statics: Melting and "f^z 

In the dilute regime, such as / < 1/24, there are numerous phase degrees of freedom per field-induced vortex 
pancake. Thus, a double transition, if there is one, might be more plausible here than at / = 1/6, one transition 
being the melting of the field induced flux lattice, the other connected to the XY-degrees of freedom p&. 

To check this possibility, we have studied / = 1/24 (a field which allows for a commensurate triangular flux lattice 
of 48 lines) on a stacked triangular grid of 24 x 24 x 24 grains. We flrst did an extensive simulated annealing run on 
a single layer, verifying that the vortices freeze into a perfect triangular lattice. We then stacked 24 such layers to 
form a three dimensional ground state. Next, the lattice was gradually warmed up in intervals of dT / J = 0.1 or 0.05, 
typically with 50,000 MC steps for each temperature. For several T close to a transition, we ran up to 10^ MC steps 
to ensure equilibration. The resulting Bragg intensity 5'(Gi) and helicity modulus component jzziT) are plotted in 
Fig. n2. [The transverse components jxxiT) and ^yy{T) fluctuate around zero for most T > 0, as expected for a very 
weakly pinned vortex lattice which is free to slide in the ah plane.] 

The results do indeed suggest the possibility of two phase transitions. The first - the melting of the vortex lattice 
- occurs near T — 1.5J = T^, where the Bragg intensity drops sharply. At higher temperatures, there is a broad 
dip in the normalized Bragg intensity which reaches a plateau at around T/ J ~ 2.1. The possible upper transition, 
near T — 2.0 J = Tg, is the point where jzziT) vanishes. Essentially the same behavior, but with an even wider 
temperature separation, has previously been observed on a cubic grid by Li and Teitel [g5| . 

We have carried out several checks to see if the separation of these two transitions is an artifact due to a finite-size 
effect. First, as shown in the inset, we monitored the dependence of jzz on accumulation time r up to 10^ MC 

steps. More precisely, we define (A\ = ^ J^ A{t)dt and use this in calculating the averages which define jzz in 

eq. (^). The < ^zz{T) >r thus defined generally evolves approximately logarithmically in r [Q until it reaches its 
apparent equilibrium value. For temperatures Tm < T < Tg, the system tends very slowly towards an apparently 
finite limiting value. We have also checked the size dependence up to 24 x 24 x 48, verifying that the 24 x 24 x 24 
behavior represents the asymptotic limit. Li and Teitel have carried out similar checks up to 200 layers in the cubic 
model |2^]. Nonetheless, ^zz has some size- dependence to a degree dependent strongly on anisotropy of the system. 

If the ratio Jz/Jxy is increased to 4.0 (where Jz and Jxy are the couplings perpendicular and parallel to the 
triangular plane), the separation {Tg — Tm)/ Jxy between the melting transition and the upper possible transition 
actually grows for a given size. For these values, the smallest-Q Bragg peak vanishes at Tm ^ 'i.QJxy, while jzz 
vanishes at Tg ~ 4.0 Jxy On the other hand, for weakly coupled layers with Jz/Jxy = 0.1, the two transitions merge 
to within less than O.lJa;^, as in the the isotropic dense / = 1/6 case. In this weakly coupled case, Tm ''^ O.GJxy 

These observations suggest that phase coherence at finite / in a disorder-free system may possibly be destroyed 
in two steps. First, coherence transverse to the average field direction is lost through melting of the lattice. But 
longitudinal coherence persists until it is destroyed, along with line-like correlations of the individual vortex segments, 
at a slightly higher temperature Tg. This is most apparent for the isotropic system only when / < 1/18. 

B. Vortex Analysis of Possible Transition at Tg 

There are several ways to look at general phase correlation function < 8(/9, z)8(0, 0) > where Q is the gauge- 
invariant local phase [gSfl of the superconducting order parameter. To probe the longitudinal phase coherence, we only 
consider c(0; z) = ^ /" P < exp[i(0(p, z) — Q{p, 0))] >, where A is the sample area - that is, the correlation function 
in the z direction. Glazman and Koshelev have pointed out that phonon-like fiuctuations in the vortex lattice lead to 
a power law decay of c(0; z) [Q (not explicitly shown in the Figure |l3|). We observe that this holds true for T < T,„. 
For higher temperatures, this dependence changes to an exponential decay c(0, z) oc exp(— z/^oz)- The correlation 
length ^oz ^ [log(^ '°^^„^//°^ )]-^ where Tq is the temperature scale such that < [6(0, z) - 9(0, z -h l)]^ >~ 1. This 

behavior is shown for T/Tm = 1-1 in Fig. |1^ for systems of several thicknesses. To ensure equilibration, we ran 86000 
MC sweeps before accumulating data over the following 30000 MC sweeps. To check the effect of boundary conditions. 



we used both open boundary (OBC) and periodic boundary conditions (PBC) along z; periodic boundary condition 
was used in the xy plane for both cases. In all cases, we observe a robust exponential dependence over a limited range 
1 < dz < £^x where deviation sets in at dz ~ ^2; ~ 12 for T/Tm = 1-1, for example. In this temperature range, we 
do not find significant dependence of c(0, z) on either system size or boundary conditions. For a given temperature 
above r„i, we can use this robust temperature regime to extract the phase correlation length ^oz fi'om our numerical 
simulation. 

The result is shown in Fig. |lj. For all T > Tm, and for separations less than £,x{T), we observe that c(0, z) ~ 
exp[— z/^zo]- £,zO gradually decreases with increasing temperature approximately as (T — Tm)~^, becoming equal to 
the unit layer spacing near T/J — 3.0 ~ Txy(O) as shown in the inset. We also observe(not shown in the Figure) that 
/ dpc{p, z) has far milder dependence on z. The exponential decay in c(0, z) is accounted for by random walk-like 
excursions of the vortex lines and the presence of dislocation or disclination loops in this temperature regime. Note 
that spinwave excitations in the vortex lattice usually lead to an algebraic decay of c(0, z). Note that given a static 
deformed vortex line configuration, we may still find a coordinate transformation {x, y, z} — > {x' , y' , z} into a curved 
space in which the vortex line is straight. In that coordinate system, we will have a long range phase coherence along 
the straight line in the z direction. Therefore, the apparent exponential decay of c(0, z) is not an equivocal indicator 
for the destruction of phase coherence along z, but gives information about the deformation of vortex lines from the 
straight configuration. 

In the large-distance tail of c(0; z), where z > l-iz/2, c(0, z) does depend on the boundary conditions and system size, 
having an upturn for the periodic boundary condition as expected. Moreover, £,zoiT) rapidly falls as T decreases toward 
TxY from above, leaving a large interval ^xiT) < dz <ti \jzI'2 in which the behavior deviates from simple exponential 
decay and is independent of the boundary condition used. We also observe that the deviations have relatively poorer 
statistics due to slow kinetics. It is likely to originate from vortex entanglement and cutting/reconnection, which 
develop on length scales larger than ^^o- Each of these rare and slow vortex crossings produces a drastic and long- 
lasting impact on the local phase correlations. The rarity of these events is due to the sizable barrier for vortex 
cutting and to the subdiffusive nature of vortex lines motion. Because of this rarity, a rapid thermal cycling across 
this sluggish region may lead to hysteresis. 

We now look more closely into the vortex configurations near Ti for / ~ 1/24. Fig. Ha shows the density-density 
correlation function n2,z(r±,L2/2) =< n^iv^^ L2/2)nz(0,0) > describing the z-component of local vorticity at separa- 
tions equal to half the total thickness Lz (=24). The most prominent feature in 71,2,2 is the disappearance of triangular 
correlations in the xy plane at melting (T^ ~ 1.5 J). This behavior is consistent with the disappearance of the Bragg 
spots in Fig. |l^. Note, however, that the central spot, corresponding to the self-correlation between the two ends of 
the same line, persists well above melting until it vanishes near T j J = 2.0, close to Tg. This is consistent with the 
fact that line-like correlations are maintained over at least 12 layers up to T/J — 2.0 as we already noted in Fig. O. 

In a sample of truly macroscopic thickness, the lines in the liquid phase should carry out random- walk-like excursions, 
leading to loss of top-to-bottom vortex density correlations over a finite correlation length denoted £^vz (T) , which may 
be of the same order of magnitude as ^oz defined above. A finite £,vz{T) means that the underlying lines are flexible^ 
not that they break up into 2D vortices. This breakup becomes relevant only for T > Tg. The objects which break 
apart into 2D objects above the melting transition are not the lines themselves, but the topological defects of the 
lattice, such as disclinations, which tend to appear as well-aligned line defects near T^ |14|j61[] . Topological defects 
look well aligned only when |ri(z) — ri{z -f l)\/aB <C 1 where ri(z) is the position of the segment of vortex line i in 
the z-th layer. Note that the relevant minimum length scale for alignment of defect lines is the mean vortex spacing, 
Ob ■ The destruction of the lattice order along the field, which may be detected by vanishing Bragg peak in neutron 
diffraction, is related to proliferation and unbinding of these defects. On the other hand, the destruction of phase 
coherence along B, as we discuss in more detail later, is related to the presence of fluctuations in the transverse 
vorticity. As the vortex density decreases, it is not a priori obvious if the energy scales for these different types of 
defects should remain the same. 

Thus the line-liquid regime, if it is really a distinct thermodynamic phase, may possibly be described as a neutral 
gas of topological defects (disclinations of both signs) within the triangular lattice in each plane, which are correlated 
over a finite length in the z direction. Above melting, one expects unbound disclinations to proliferate. Hence, the 
long-range structural correlations of the vortex lattice are lost in all directions upon melting. However, the phase 
rigidity, as measured by 7^2, may persist even above melting, but scaled down by the factor of ^^o/Lz, the fraction of 
the volume of the sample into which the applied twist penetrates. Presumably, this continuous suppression, unless pre- 
empted by a first order transition (for high densities), persists until the condition Lz/^zq — > 00 is met via proliferation 
of "unbound vortex loops" (vortex lines extending an infinite distance in the transverse direction). At this point, the 
phase coherence even between neighbouring planes normal to B will be lost. 

In Fig. nq, we show snapshots of vortex configurations at T^, Tf and a temperature between T^ and Tg. In this 
regime, by using a bond-searching algorithm, we have identified three distinct classes of vortex lines. The first consists 
of small vortex loops which close on themselves without crossing either of the two opposite bounding surfaces. The 



second class contains all isolated lines beginning at the bottom xy plane and ending at the top one. Most of the 
disentangled field-induced vortex lines fall into this group. Finally, there occur "vortex tangles". These are lines 
connected at a given time to one another by the crossing of two vortex segments in the same unit cell. Such tangles 
are formed either by collision of two flux lines or by interactions of such lines with the vortex loop excitations. This 
tangle is not static: the collisions which produce it are more and more frequent with increasing temperature and its 
overall shape will evolve with more rapidity as T increases. 

The three columns of the Figure represent the fraction of the vortices belonging to each class at a given instant. 
On melting {T / J = 1.5), the fluctuating lines in our finite sample still remain largely disentangled and separated 
from each other. As T ^- Ti, the density of loop excitations increases (left column), while the field-induced lines 
(central column) have stronger lateral fluctuations . Both of these effects cause more and more "connected" clusters 
(i. e., vortex tangles) to appear. Finally, at T^, an infinite tangle, connected by crossing vortex lines, forms. At this 
temperature, the connected tangle of vortices form a (D — 1) dimensional manifold of a tortuous shape, transverse to 
B, and cut the original D{= 3)- dimensional coherent XY system into halves. 



Fig. |17| shows an instantaneous vortex cluster size distribution for various temperatures at / = 1/24. To generate 
this distribution, we define the projected transverse length of each vortex loop (or tangle) by Ixy = ^ |z x ni,| for each 
isolated cluster composed of unit vortex segments rii,, and accumulate a histogram, V{ixy)- [We consider only the size 
distribution projected onto the xy plane because the field induced lines (which are infinite along z) could mask the 
loops with large extent in the z direction. We also believe that these fiuctuations are more relevant to the vanishing 

Oi'fzz-] 

In the first panel of Fig. n% we plot V{ixy) for several T < T,,,,. Each plot has a sharp maximum cutoff and 
a pronounced peak, which is due to the finite average lateral fluctuations of the field-induced vortex lines. For 
r,„ < T < Tg (second panel), the weight of distribution is shifted towards larger i^y, because lines in the liquid phase 
undergo larger transverse fiuctuations. Closed vortex loops also begin to appear in this region. As T increases, the 
distribution is cut off at progressively larger values, as more and more lines join the connected clusters. Finally, for 
T > Ti, all curves are characterized by the appearance of "infinite" clusters, with no obvious length cutoff. The 
distribution appears to fall off algebraically in this regime - i. e., Viixy) ^ ixy with fi ^^ 1.0 < 2 - suggesting 
< ixy >^ oo for T > Ti. 

In Fig. |l8|, we show the maximum value ixy occurring over 10, 000 MC sweeps for each temperature. The size is 
normalized by the linear system dimension in the xy plane (48), and also by the number of xy- planes (24). Although 
£xy grows monotonically for T > T„j, it seems to jump discontinuously from ^ 1 to ^ 2 between T/J — 2.0 and 2.1 
(near Ti for samples of this size). Qualitatively similar behavior occurs for the isotropic f — XY transition near 
TxY{f = 0) [cf. Fig.|. 

C. Entanglement, Winding Number and Other Exotica 

We now discuss a possible extension of the vortex loop picture of the zero field XY transition to the hypothetical 
Tg transition or crossover for / < 1/24. Such loop excitations have received far less attention in 3D systems pi] than 
their 2D counterparts, possibly because they require more energy to excite and therefore matter only very close to the 
mean field transition. But in high-Tc materials, the short correlation lengths, high anisotropy, and high Tc broadens 
the vortex- loop-dominated regime ||62| , |63| , before amplitude fluctuations set in. 

In a cubic sample with periodic boundary conditions, all vortex lines naturally close on themselves to form loops. 
These loops are of two topological types: those which can continuously shrink to a point ("trivial class") and those 
which cannot ( "nontrivial" ) . The latter are said to have a nonzero "winding number," i. e., number of infinite lines 
in a given direction. In the 2D periodic case, the loops lies on the surface of a torus. Here, there are two distinct 
subclasses of non-trivial loops: one which winds around its circumference, and another which runs transverse to it. In 
the infinite 2D geometry, these correspond to lines infinite in cither the x or the y direction. On the 3D hypertorus, 
there are infinite lines in any of three directions. 



These notions play critical role in the description of dissipation via vortex motion, i. e. phase slips |64|. For current 
flowing in a given direction, the dissipation may occur either through the expansion of loops, or through motion of an 
infinite line. In either case, the dissipation arises from fiuctuations in the winding number of vortex lines perpendicular 
to the current(c. f. Equations (||) and (jj)). Note that when a finite field is applied along z direction, the hypertorus 
already contains many "windings" along that direction even without an applied current. 

In the absence of pinning, dissipation in the plane normal to z is governed by fiuctuations of the winding number 
in the z direction. This dissipation should not depend directly on whether or not the "windings" - that is, the vortex 
lines - form a lattice, but may depend on the density and mobility of windings. 



Dissipation parallel to the field direction (c-axis resistance) depends mainly on winding number fluctuations trans- 
verse to z. Clearly, the average winding in this direction vanishes, unless the field-induced lines themselves, while 
winding along z as required, also wind along another direction like wires around a solenoid. For this to occur, the lines 
would have to break a chiral symmetry, spontaneously generating a global surface current with a net magnetization 
normal to z - an effect which should be prohibited energetically in the ground state. 

It may occur, however, if there exist entangled field induced vortex lines which collide with each other to switch 
connections(a process we may call "cutting and reconnection" ) . In Fig. n9, we show two field induced lines residing on 
the surface of a torus (left panel) going through such a cutting and reconnection [(panels (a)-(b)]. The right column of 
the figure shows an alternative view of the same process in an infinite space with open boundary conditions. Initially, 
both vortex lines wind only along the z axis. After the cutting and a special reconnection process in which one 
strand circles around the torus before meeting its other end, a net transverse winding number has been created. This 
"global" process is, however, energetically expensive because it involves a spatially extended excursion [panel (b)] and 
should occur very rarely, even in the melt. 

If we introduce an "entanglement length" ^c, defined as the average distance along z required for any two vortex 
lines to wind around each other, we expect tc to be infinite for T < Tm, but to become finite in the liquid phase. 
Because of the finite line tension and repulsive interactions between vortex line segments, such entanglement events 
along the flux lines are costly in energy and hence rare, in the liquid near melting. Deeper into the liquid phase, as 
the repulsive interaction between vortex lines is overcome by entropic forces of attraction, £c should become much 
shorter, leading to a much denser entanglement pattern. The now numerous local transverse fluctuations, and local 
cutting and reconnection events (i. e., collisions) generate fluctuations in the "global" transverse winding number and 
cause jzz = 0. On symmetry grounds, the average transverse vorticities < n^ > and < Uy > have to be zero at all 
temperatures. But possibly < \nx\^ -\- \ny\^ > acquires a finite value for T > Tg, suggesting that this quantity could 
be used as another "order parameter" for the hypothetical phase transition at T = Ti (with a nonzero value at higher 
temperatures) . 

Alternatively, we may view the upper transition in the context of a bond percolation transition. The field induced 
lines provide a kind of backbone network. With increasing T, vortex lines undergo more and more transverse collisions. 
At T ^ Tg, these collisions induce the entire ensemble of field-induced vortex lines to form an infinite connected D — 1 
dimensional structure transverse to the applied field, causing large fiuctuations in the transverse winding number 
[thick gray line in panel (c)], thereby wiping out any superconducting path connecting the top and the bottom layers 
normal to the field. 

Let the mean-square transverse displacement of field-induced vortices per layer be denoted £^ =< \ri{z) — ri{z — 
1)1^ > . Then ic/d is defined as the number of layers along z over which a line wanders transversely by the average 
intervortex distance. We write this condition as £|n • [Ic/d]'^'' = a\, where we introduce an unspecified "wandering" 
exponent C,. In the limit of dilute (independent) lines, we expect C, ~ 1/2, corresponding to a random walk of each 
vortex line segment. Long-range intervortex repulsion is known to renormalize the unit step size £t from c( T)-{T / J^)^'^ 
down to a smaller value with a similar form with an unspecified T dependence encoded in c{T) < 1 [£^,^. It is 
not clear how C is affected by intervortex repulsion, but possibly the interactions with other fiuctuating lines are 
equivalent to the line of interest being in a random environment. For a fiexible line in a 3D random environment, 
C "^ 0.6 [^^lial- ^'^^ D > 2, in the actual system of interacting fiuctuating lines, no exact result is available for C, 
although several numerical results and conjectures give C, ~ 0.2 — 0.6 169]. 

We can use these crude estimates to make a guess at the field dependence of T^, interpreted as a bond percolation 
transition. Along a given field-induced vortex line, the probability per unit length that a transverse connection is 
made to a neighbouring vortex line at any position along the z-axis is p = d/lc = (^r/as)^^''- Since €t oc {T/Jz)^.5 
and qb oc B~'^-^, the percolation threshold is reached roughly when c{T)'^TB/Jz > [Pc]^^, where Pc is an appropriate 
percolation threshold. This condition defines a lower bounds for a possible transition at Bg{T) which approximately 
follows 

For a dense lattice or large anisotropy (i. e., small Jz), this condition is probably satisfied immediately upon melting, 
as at / = 1/6. For dilute systems, however, the second transition is not automatically triggered by melting and may 
occur only deep into the liquid phase, at a temperature where the entanglement barrier is sufficiently weak to allow an 
infinite vortex tangle to form. Whether this percolation transition is a true phase transition or only a sharp crossover 
remains to be determined. 



This same picture hints at how correlated pins such as columnar damage tracks |70 may increase Ti. Such columnar 
disorder will encourage the vortex lines to stay straight along the defect track, reducing the effective unit step £t by 
a factor Cp <C c. As a result, the wandering exponent ^p may also change from its thermal value ^. Consequently, Ti 
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will be enhanced by an overall factor of ( — )^(— )^^^ ''''•'• 

In summary, the upper transition is characterized by the following set of equivalent criteria: pq , p3| Disappearance 
of finite transverse diamagnetism; disappearance of phase rigidity along the field direction; appearance of an infinite 
transverse vortex cluster; large fluctuations in the global transverse winding number or the net vorticity A4; onset of 
finite c-axis phase-slip resistance in the limit of vanishing bias current in the c direction, which is equivalent to saying 
no superconducting path exists over macroscopic distances. 

D. Dissipation for / = 1/24 

While Bitter decoration serves as a detailed probe of spatial vortex configurations ||6l|,|7l|,^ , it yields ambiguous 
information about freezing, and is restricted to very low flux densities. Cubitt et al obtained evidence of a melting 
transition in Bi2SrCa2Cu208 from low angle neutron diffraction 0]. Similar results were obtained by a fiSK technique, 
which probes the local magnetic field distribution g] . NMR [Q and atomic beam [Q techniques have also been used 
to study both the static properties and the melting of the vortex lattice. More recently. Schilling et al employed a 
differential thermometry to search for the latent heat of melting in YBa2Cu307_5 and to obtain a melting curve [^,0. 

Far more information has been accumulated from transport measurements, but this is much less easily interpreted 
in terms of vortex lattice melting. The interpretation is complicated by disorder, as well as by the fact that the 
measurements are nonequilibrium and usually involve nonuniform current distributions. Safar et al [pj measured a 
sharp jump in resistivity in the mixed state of YBa2Cu307_5 . The resistivity also showed a hysteretic behavior upon 
thermal cycling, indicating a first order transition. The transition line thus obtained seems to coincide with "melting 
curves" obtained by torque measurements [^ , and more recently, by differential specific heat measurements g] . Kwok 
et al have carefully demonstrated the effect of twin boundary pinning on the melting transition in a series of transport 
measurements which track the so-called "peak effect" associated with vortex lattice softening |2^| . They find that the 
peak effect sets in at a few degrees below the melting curve determined from a sharp kink in resistivity. This sharp 
resistivity kink, as observed by both Safar et al B and Kwok et al |7q|, tends to become less pronounced both at very 
high(B > lOtesla) or low fiux densities(i3 < Itesla) @. 

An ideal, but impractical, transport measurement to determine the melting curve would consist of applying an 
infinitesimal current to induce a net Lorentz force on the lattice, which is held in place by a balancing pinning force. 
As soon as the lattice melts, individual lines would begin to drift, inducing "flux flow" resistance. Most real materials, 
however, are complicated by disorder, and even the static properties of the lattice with disorder are incompletely 
understood f7^,Q. In the presence of disorder, varying the field density produces changes in both the effective 
pinning strength and the effective flux lattice anisotropy. Depending on relative strengths of all these cornpeting 
effects, many complications may arise in probing thermodynamic properties using transport experiments [[79|-p3|. 

A number of recent transport measurements have, nonetheless, produced rich information about flux lattice melting. 
Pastoriza et al have used a non-uniform distribution of pinning strength to probe the shear modulus directly [[44| , 
giving direct information about the lattice stability. Zeldov et al [gj have used local Hall probes in Bi2SrCa2Cu208 to 
monitor the local field density. They found a very sharp jump, again interpreted as a signature of a first order 
melting transition. More recently p3, Fuchs et al performed simultaneous measurements of the resistance and 
local magnetization, confirming that that the jump in local magnetic density coincides with a sharp increase in 
resistance (but not necessarily a discontinuous jump) . The heat of melting per vortex per layer inferred from this local 
magnetization jump, however, shows rather peculiar features: it vanishes continuously as the field is increased, while 
steeply increasing as the field is lowered toward zero. 

These results raise several outstanding questions: How can the seemingly first order transition line terminate 
apparently at a point in the H-T plane? Does the melting line monotonically approach Tc{H = 0) or follow a 
reentrant melting curve? Another important issue is the longitudinal phase coherence probed by c-axis resistivity 
vs. T measurements [p4| , p5|j65| ] , which show a striking series of broad peaks in Bi2SrCa2Cu2 08 single crystals. The 
nonlocal conductivity associated with this phase coherence can be probed in the so-called flux transformer geometry. 
The experimental data of Keener et al [Q showed that phase coherence over a flnite correlation length along B 
persists above the melting transition of the vortex lattice in some region of the H-T phase diagram. 

A simple and natural model for probing the dynamics of the mixed state is a network of resistively shunted Josephson 
junctions with Langevin noise. In this section, we present some results of simulations using this model, and to connect 
these to the analogous static XY results. Our calculations are carried out as follows. At any given temperature, the 
final snapshots from the MC simulations are used as the initial dynamical phase configurations. We use an integration 
time step At = Q.Uq. After the current is switched on, 1000 — 5000 Ai is allowed for the system to reach a steady 
state, following which the voltage is averaged over the next 6000 — 12000 steps of Af . 
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Fig. po| shows the "bulk in-plane resistance" at / = 1/24. The measurement geometry is as shown in Fig. |] (a); thus, 
these calculations probe the shear rigidity of the lattice in contrast to the usual transport experiment in which random 
pins play an essential and complicating role. Through the rest of this paper, we will call these calculated quantities 
Rab and R^ even though they differ from the resistivity measured in most transport experiments. At the highest bias 
current of 2.83/c per grain(equivalent to 1.4/c per bond), we have a smooth curve without any noticeable changes 
either at T„j or at Tg. For lower values of driving current, sharper features emerge. There is a slope discontinuity 
near Tm = 1.5 J for both I/Ic = 0.83 and 0.083, but the most dramatic change occurs near Tg = 2.0 J. The entangled 
line liquid for T^ < T < Ti seems to have a sizable viscosity. This viscosity impedes the motion of the flux lines in 
the liquid, the two halves of which are driven past each other by opposing Lorentz forces. As a result, the lines move 
slowly, and dissipation (defined as a "resistance" Rab) is small. The steady increase of Rab with temperature in this 
region is due to screening by vortex loops which gradually lowers the viscosity. For T > Ti, this viscosity vanishes, 
leading to a steep increase in Rab This increase at Tg is enhanced by additional (and probably dominant) dissipation 
produced as the system goes through an XY-like transition or crossover. The large viscosity for Tm < T < Tg is also 
consistent with the slow (In t) equilibration seen in the Monte Carlo measurement of jzz for Tm < T < Tf [Q . We 
believe that the change in Rab near T( shares the same mechanism as that seen near Txvif — 0) shown in Fig. g. 

Fig. pl| shows the "c axis resistivity" i?c at f = 1/24, as calculated using the geometry of Fig. n^ (b). For comparison, 
we also show the calculated ^zz- As the driving current is reduced. Re seems to approach a curve which vanishes 
asymptotically as T ^ T^ , coinciding with the vanishing ^zz- Our numerical results thus suggest that the dramatic 
increase in R^ results from the Tg transition rather than melting. The increase in R^ is thus correlated with massive 



vortex line cutting, as put forward earlier |p6 65 1, and with an increase in the density of transverse vortex segments 
< \nxy\ >, the frequency of vortex line crossings, and fluctuations in the transverse net vorticity, SAd^. This distinction 
between T^, and Tm may be most important for Bi2SrCa2Cu208 at low fields, and in disordered dense systems [ p3[ , 
where the temperatures may be most separated. 

We can draw some conclusions relevant to experiment from the current dependence of both Rab and Re observed in 
our simulations. First, the "melting line," as detected via a voltage criterion at constant current in a bulk resistance 
measurement, should be sensitive to the driving current, even at a very low bias. Existence of pinning force is essential 
in getting distinct transport behaviors for the lattice and the liquid. On the other hand, the transition at Ti, whether 
monitored by the vanishing of Re in the limit of small current or by a jump in Rab between two finite values, should 
be relatively insensitive to applied current density, since the main mechanism of dissipation(presence of transverse 
vorticity) in this case is switched off below Tg and sets in above Tg irrespective of whether we have pins or not. Indeed, 
just such an observation has been made by Keener et al [B2| in describing their curves for Tm{H) (melting) and 
Td{H) ("decoupling transition"), as obtained by flux-transformer measurements on Bi2SrCa2Cu208 single crystals. 
It is plausible that their Td{H) at very low fields {B < 100 Gauss) corresponds to Tg in our model. True melting 
line is presumably the limiting value of the current-dependent Tm{H, J) as J ^ 0. It is not experimentally verified 
whether such a limiting value coincides with Td or not. 

Let us briefly comment on the experimental possibility of distinguishing between Tm and Tg If we imagine a 
hypothetical isotropic high temperature superconductor, the coupling constant J in our model calculation is related 

to the parameters of the superconductor via J ~ igTr^A" (o) ^ (-'^ ~ [^/^co]^) (assuming the two-fluid model). Taking 

Tco = Q2K, d = lOA, and A(0) = lOOOA, we find Tm ~ 89.7K and Tg ~ 90.3ii: [eq. (11)]. Thus the two transitions 
are remarkably close even for the isotropic case. In real materials such as YBa2Cu307_5 and Bi2SrCa2Cu208 , the 
separation between the two will be further reduced by an anisotropy factor, although pinning disorder may tend to 
separate them. Therefore, in many cases, it will be nearly impossible to separate the two transitions experimentally. 

VII. LOCAL MAGNETIZATION JUMP AND HEAT OF MELTING 

A striking result of the Bi2SrCa2Cu2 08 micro-Hall probe measurements is the sharp jump in local magnetization 
(vortex density) across the phase transition |0|. At "high" fields (~ 2000), the jump occurs at nearly constant T, 
and even for lower fields, still within 5T ~ 3mK. The heat of melting per vortex per layer, Tm/^S = — '-^ — ^^j 
as obtained from the Clausius-Clapeyron relation, increases monotonically from at -B ~ 4006* to about O.Gfc^ at 
B ^ 55G, beyond which the slope ^^P^increases very sharply (cf. Fig. 6 of g]). 

Similar jumps also seem to occur in YBa2Cu307_5 , as reported in recent calorimetric B and magnetization 
measurements ^ , p7| . The estimated latent heat of melting yields AS* (per vortex per layer) ^ OAks for 1 < B < 
8(tesla). The data (cf. Fig. 1 of Ref. @) shows that the jump AM at T = 85K {B - 3.7 tesla) for YBa2Cu307-A- is 
spread over a field range 6B ^ O.IT, or, for a given field, over a temperature range ST ~ O.IK. This jump decreases 
rather abruptly for flux densities B < IT. The estimated entropy of melting (~ O.AkB per vortex pancake) is quite 
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close to that numerically obtained by Hetzelei al ^^ (~ O.Sfcs per pancake), and also to the values obtained in model 
calculations based on the lowest Landau level and London approximations |88[| . 

As the field decreases, the jump in M occurs over a broader temperature range (cf. Fig. 3 of |^). The resistance 
jumps measured by Kwok et al |75| ], attributed to the melting transition, also become broader with decreasing field. By 
contrast, the height of the resistivity kink seems quite uniform over a wide range of fields. These last two observations 
are consistent, however, if we interpret the jump in resistance as a signal that 7^2 —f 0. In view of all these facts, 
it is plausible that, at least at relatively high fields which corresponds to fJxy/Jz > 1/18, the experimental jumps 
observed in local magnetization and resistance |3|,|6|,|75| shows the combined effects of two distinct processes, occurring 
within AT < lOmK. As a corollary, the very low field measurements(/ J^y/Jz < 1/18) actually may not track the 
melting transition itself, but various manifestations of the predominant XY fluctuations (vortex loops) which are most 
conspicuous near T^.. 

Observation of the reentrant melting curve has not been reported in any high-Tc materials. In NbSe2, the melting 
line detected by the peak effect was reported to be non-monotonic in field [pO|, consistent with the reentrant melting 
curve proposed for the more anisotropic high-Tc superconductors [^. But if one interprets the magnetization jump 
in Bi2SrCa2Cu208 as evidence for melting, then the melting curve for Bi2SrCa2Cu2 08 apparently approaches Tco 
monotonically at field as low as ~ IG. This behavior is surprising since, at these fields, the vortex separation far 
exceeds the magnetic screening length. Furthermore, in YBa2Cu307_5 , the peak effect at 0.35 — 1.5 T is observed 
to lie below the resistivity kinks [£1| (about G.SK below at 0.5 T). These measurements suggest that, at least for low 
flux densities, transport measurements may actually not be probing flux lattice melting. 

We propose that low-fleld melting is indeed reentrant. Most low-fleld experiments which probe magnetization p,p4[, 
thermal properties [||, and transport coefficients |8| p^j3l| , [75[ actually track T^, which is progressively more separated 
from Tto and approaches Txf(O) as the held is reduced. We have already shown that dramatic changes in Rab and 
Re occur at Ti. Moreover, the broad peaks in Cy are centered at Ti and they, too, approach T'xy(O) as B decreases. 
Our estimated upper bound for the total entropy release in the temperature range T^ < T < Ti, as estimated from 
the T— dependence of the internal energy, qualitatively resembles that of Zeldov et al (Fig. 6 of 0) in that it steeply 
increases as / decreases. At higher flelds, the observed (and also calculated) AS* ~ 0.3 — O.bks is consistent with the 
destruction of phase coherence at a single transition. By contrast, at low fields, phase rigidity is lost in a two-step 
process. Most of the entropy release {AS > O.^ks per vortex per layer) occurs near T ^ Ti> T,„, whether or not this 
is a true phase transition. Note that a hysteresis in the resistivity may be observed near T£ due to flnite vortex-cutting 
barriers below Tg. This is not necessarily an evidence for a flrst-order melting transition at very low fields. 

The melting line in the dilute limit may be quite difficult to detect experimentally. Conceivably it may be tracked 
by the peak effect, by high-resolution IV measurements, or by direct measurement of the shear modulus [Q. Of 
course, direct observation of a vanishing neutron diffraction pattern as in H] would be ideal, but this technique is of 
limited applicability in this density range. 

To shed further light on this problem, we have carried out calculations with mixed boundary conditions. That is, we 
allow local density fluctuations in the net z-component of vorticity by using free boundary conditions in x~ and y— 
directions, while retaining periodic boundary conditions along the z— axis. Of course, surface effect are now stronger, 
possibly reducing the melting temperature. Another point is that our uniform- frustration model assumes that A = 00. 
Therefore, we should proceed with some caution in relating our numerical results to experimental data. 

To study the system with these mixed boundary conditions, we again did a simulated annealing run for a single 
layer of the triangular grid to find the lowest energy configuration. By stacking the resulting state layer by layer, 
we form the ground state lattice, which, because of incommensurability and the free boundaries, now consists of an 
imperfect triangular lattice with some defects. This lattice melts at T/ J < 1.4 for the nominal density of / = 1/24 
on a 26 X 26 X 12 grid. This is slightly below the value Tm/ J ~ 1-5 found for the fixed density system of 24 layers 
with periodic boundary conditions. For a nominal density of / — 1/6, melting occurs near T/ J < 1.15 with these 
mixed boundary conditions. 

One might think of defining the "magnetization" M^ as the average net vortex density n = J nz{r)dr/A, where 
nz{r) is the local vortex density and A is the total area. However, M^ defined in this way suffers from spurious 
boundary effects, arising from the depletion of vortices near the boundaries in the lattice phase |p9]. Upon freezing, 
the lattice develops a rigid iree surface of irregular shape, expelling some of the vortices from the rectangular bounding 
box. The resulting change in density, Sn/n, is an artifact of the open boundary conditions, and we find that it vanishes 
for large samples as 1/vA, confirming that it originates from a surface effect. 

Instead, we define Mz by a criterion involving the local Voronoi cell area At, i. e., the area of the generalized 
Wigner-Seitz cell for vortex i (the shaded area shown in Fig. k3). Before applying the procedure, we first eliminate 
the thermally induced vortex loops, which are present in addition to the field-induced vortices for T > 0.5 x Tg. To 
do this, we pair each antivortex with the nearest vortex in each plane, identifying the resulting pairs as bound dipoles 
to be excluded from the count (cf. left panel of Fig. E3h. Since most such dipole pairs have linear dimensions much 
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smaller than \/\/< n >, this criterion is justified. We then perform a Delaunay triangulation on the field induced 
vortices to determine topological neighbors for each vortex. From the bond configuration thus determined, we obtain 
its dual, which is the desired Voronoi diagram. A local vortex density at a point R may then be defined as 

n(R)=^(5(Re A)/A (10) 

i 

where (5(R G Ai) = 1 if the point R lies in the Voronoi cell associated with vortex i, and zero otherwise. Next, the 
local magnetization M^, which we interpret as the hulk average density < ti > is calculated from 

<">=Ex' (11) 

i. e., as the average of the inverse Voronoi area for vortices lying within a measurement area C suitably distant from 
the sample boundary. 

In Fig. gj, we show the relative average vortex density (filled circles) along z, < n^ > /tiq, normalized to the 
nominal density per layer at / = 1/6 and / = 1/24. For / = 1/6, < n^ > shows a sharp jump at T„ to a value about 
7 % larger than tiq. For / = 1/24, the there is a similar change of about 15 % which is less sharp than at / = 1/6 
and is centered at Tg. From these two data points, we observe that [< n > (T^) — no]/no ~ /~^^^. Comparing the 
result for two different sample areas L^jLj, for / = 1/24, we have verified that the observed change is indeed a hulk 
phenomenon, independent of any surface influence. Note that we have about the same number of flux lines(~ 0(200)) 
for both / = 1/6 and / = 1/24 for our chosen sample sizes of 24 x 24 x 12 and 48 x 48 x 12. While the jump 
occurs at Tm for / = 1/6, we do not observe a similar feature near melting {Tm = 1.35J) a,t f — 1/24. Therefore, the 
cause of the jump in the local vortex density should be sought in the nature of transition at Ti, rather than in the 
mechanism for flux lattice melting. Note that these jumps resemble those in Fig. 5 of Q in the "anomalous low-field 
regime" {1 < B < 55G) in the following sense: the fractional change in vortex density decreases, and the jump 
becomes sharper, as the field increases. Probably, the line density (Hz) increases with increasing T for T^ <T <Ti 
because the repulsive intervortex interaction is screened by polarizable vortex loops. The 2D analog of this effect is 
the screening of the repulsion between field- induced vortices by thermally excited vortex-antivortex pairs pffl . 

Does the jump in flux density occur exactly at the melting transition, or is it more closely connected to the other 
"transition" at Ti, i. e., to a transition between two liquids with different compressibilities? This question may 
actually be rather academic, since Tm and Ti may practically merge in real, anisotropic materials at high fields. Eq. 
(0) gives a rough criterion for Ti in isotropic systems: 1/lc = {c^BT / JzY/'^ > pc. For anisotropic systems such 
as Bi2SrCa2Cu208 and YBa2Cu307_5 , a given value of / corresponds to a field which is reduced, relative to the 
isotropic system, by a factor of Jz/Jxy Therefore, the merging of Ti and Tm, which in isotropic systems occurs 
around / ~ 1/18, should in anisotropic materials occur around f = (1/18) J^/J^y, This anisotropy factor Jz/Jxy could 
be as small as 0(0.0001) in Bi2SrCa2Cu208 . 

VIII. DISCUSSION 

A. Analogy to XY Transitions of Slabs of Finite Thickness 

In the previous sections, we made following observations from numerical simulations: (i) ^zz vanishes at Ti ^ Tm 
for isotropic system with / > 1/18 and Tt{B) appears to terminate at Txy for B — > 0; (ii) As B changes, it tracks 
the broad peak in Cy which behaves 

cr"(i?)~-MLB), (12) 

and 

|Tf-rxy|^L-'/', (13) 

where Lb = /~^^^ ~ B^^l"^; (iii), Ti, and not T^,, appears to coincide with the principal change in local vortex 
density which follows 

AA//B - B^^l'^. (14) 

There is a corresponding change in bulk resistances Rab and Re over the same temperature range; (iv) T^ and Tm 
seem to merge at a sufficiently high field. 
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We now summarize som.e recent experimental observations which appear to be consistent with these numerical re- 
sults. Schilling et al |7[] have reported high resolution calorimetric evidence for a first order transition in YBa2Cu307_5 , 
which they interpret as a melting transition. They observe a very sharp delta- function like peak lying on the left 
shoulder of the broad peak in specific heat in the range of 0.75 — 9 tesla, which roughly corresponds to 1/81 < / < 1/6 
in our isotropic sample. The delta-function appears to vanish for densities lower than about 0.5 tesla. This re- 
markable experiment thus establishes the existence of a first order melting transition line, which empirically follows 
\Tm — T'c(0)| ~ Lg^'^^ over 0.75 — 9 tesla. These data are consistent with our numerical results on the following points: 
(i) a first order melting transition exists, and becomes weaker as the fiux density is lowered; (ii) the melting transition 
is located on the left shoulder of a broader peak in Cv ', and (iii) the height of this broad peak and its position generally 
follow the behavior described in equations |lj, |l^ and |l^. As further evidence of the correspondence, we show in Fig. 
p5| , the data extracted from Fig. 1 of H] for fields as high as 6 T. At higher fields (> 7 T), the points deviate from 
the observed power law behavior shown in the Figure. 



Another experimental data which agrees well with our numerical results is that of Roulin et al 1 59 1 . These workers 



have reported that both the melting curve Tm{H) and a point they label the "superconducting- normal" (SN) transition 
as monitored by tracking the maximum in Cv as a function of H both follow the equation \T — Tc(0)| ~ Lg , 
consistent with our numerical results and the scaling analysis discussed above (to within logarithmic corrections). 
Welp et al [p^ have reported a detailed study of AM as a function of T and H. The data presented in Fig. 2 of 
their paper shows that AM/B ^ B"^/^ for 1.8 < B < 5.6 tesla, once again in agreement with both our numerical 
results and the scaling data From this data, we conclude that our numerical observations, based on a frustrated 3D 
XY model, are generally consistent with recent experimental observations. 

Most interpretations of these experimental results focused on only one true phase transition in the low field regime, 
namely a first-order liquid-solid transition. This viewpoint is consistent with our numerical results only if we assume 
that Te — Tm (where Te is defined as the temperature where jzz vanishes) will go to zero in the thermodynamic limit. 
In the following, we will briefly review a multicritical scaling approach which assumes a single melting transition which 
happens to be in the vicinity of the zero field XY critical point. Our data are not sufficient to determine without 
ambiguity whether or not this assumption is correct. Therefore, we follow it by giving an alternative discussion based 
on the hypothesis that there are actually two separate phase transition lines: Tm{H) for flux lattice melting and 
Ti{H) for complete destruction of any superconducting path (phase coherence) in all directions. 

Friesen and Muzikar [Q describe the SN transition at a finite B in the vicinity of the / = XY critical point. 
Their scaling hypothesis takes the form 

/.(i3,T)^itr"0±(i?ir''') (15) 

for the singular portion of the free energy density in the XY critical region. Here a and u are the standard critical 
exponents describing the specific heat and correlation length of the / = critical point, t — T — Txv'(O), and (j)± 
are appropriate scaling functions. B is put in by hand based on the assumption that it is the only relevant length 
scale. It is plausible, but does not have rigorous justification. Since a ~ and v ^ 2/3 for the d = 3 XY model, this 
expression can be rewritten as 

/.(B,T)^|i|2ln|i|0±(i3|r'/')- (16) 

The singular part of Cy ^ —d^fs/dt^ can now be shown to satisfy the relation (for T < Txy) 

Cv{B,T)r^C{x)lnt (17) 

where x = i?|i|~^/^ is the appropriate scaling variable and C{x) is another scaling function. From this we find (i) that 
the quantity Cy/ In |i| has a maximum at some fixed value of a;; and (ii) at that fixed value of x, the maximum value 
(jrnax ^ \ii^t\. Both (i) and (ii) are in agreement with our numerical data. Similarly, the magnetization is given by 
M{B,T) - {df/dB)T. It is readily shown to satisfy 

M r-. M{x){B^I^\n\x\-\nB), (18) 

where AA is another scaling function. A reasonable interpretation of the "jump" AM in magnetization is the difference 
in M between two fixed values xi and X2 of the scaling variable. Then, if the term involving In B can be neglected, 
we have l^M / B « B~^'^ in agreement with our numerical results. 

This same scaling picture can be used to interpret the heat of fusion at the first-order melting transition at Tm{B). 
Assuming that T„i{B) happens to be in the XY critical region, we write the free energy densities below and above 
Tm{B) as -t^\-a.\t\fXB\t\-'^'^) and -t^\a\t\fi{B\t\-'^/'^), where fs and fi are two different scaling forms for the 
free energy density above and below the melting transition. At the melting point, these free energy densities must 
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be equal. Then a little algebra shows that the jump As in the entropy density S = —{df/dT)B takes the form 
As = —t^ln\t\{f'^{xm) — fi{xm)), where Xm = B\tm\~^^^ is the value of the scahng parameter at the melting point. 
Using B ~ \t\^^^ along the melting curve, we find that along the melting curve 

Asr^B^/^lnB. (19) 

Thus the melting transition should have an entropy jump which gets smaller as the field is reduced. 

If T£{B) represents a true phase transition as we desribcd using the idea of vortex tangles, how can it be understood 
in terms of the phase coherence? One possibility is a line of critical points for a continuous phase transition similar to 
the XY transition in a semi-infinite slab. This view provides a natural explanation why the scaling theory with the 
scaling variable B£^^ should be successful. Mathematically, one can attach a "phantom" cut-line to the core of each 
vortex segment across which phase slips by 27r. These are benign since continuity and single valuedness of the phase 
& at every point is ensured. However, their shape and motion can be monitored most conveniently to keep track 
of the spatial and temporal disturbance of the phase coherence which have important consequences such as phase 
slip dissipation in superconductors. For an isolated vortex segment placed at origin with positive vorticity along z, 
the cut-line may lie straight along the positive £-axis. A negative vortex will then have the cut-line on the negative 
X-axis. Once we choose a cut-line for a particular vortex by fixing the reference phase angle, it provides the reference 
for all other vortices. Note that these cut-lines can only terminate either at the sample boundary or at the core of 
vortices of opposite charges. When there are several interacting vortices, these cut-lines are no longer straight, and 
their tortuosity reflects the phase disturbances induced by deformation of the vortex configuration away from perfect 
lattice. If the vortex lines were straight, we will observe that the cut-lines associated with each segment all line up as 
we move along a vortex line. Therefore, a cut-line associated with each vortex line will form a semi-infinite cut-sheet, 
separated from other sheets by roughly Lb ~ B~^/^. As the vortex lines become tortuous in the liquid phase above 
Tm, the cut-sheets will become wrinkled and our system will look like a three dimensional maze walled by these 
sheets. Both in the vortex solid and line-liquid phases, this maze will allow a arbitrarily curved path connecting both 
sides (either along z or i-axes) of the sample, and the average width of the path free of the walls will be 0{Lb)- We 
conjecture that it is possible that the phase variables may maintain long range coherence. At low fields(large Lb), 
this tortuous slab contains many XY phase degrees of freedom which, being confined within the walls, do not feel 
the presence of free vortices, and therefore could conceivably undergo a phase transition in the universality class of a 
zero-field XY model in a "film" of thickness ~ Lb, i. e., a quasi 3D-XY transition. This crosses over to a bulk XY 
transition as i? ^ 0. 

There are two possible objections to this picture. First, our numerical results only hint at, and certainly do not 
prove, two separate phase transitions. Secondly, the "film" mentioned above is a dynamical rather than an equilibrium 
film, in the sense that its boundaries are not fixed. It is not clear that such a dynamical object could have an XY 
phase transition. The boundaries (i. e., cut-sheets of the deformed vortex lines) are, of course, moving subdiffusively 
|0,p2| as long as S,vz/d > 0(10). This condition, as we confirmed numerically in section VI B, holds true in the range 



Tm < T < Tg and makes the above picture more plausible. It also greatly enhances the chance if we consider pins in 
real material, since even a single vortex line then becomes collectively pinned into a glassy state. 

We now discuss a 3D XY-like transition for the infinite slab of thickness Lb- Such a slab belongs to the G'2-class 
in Barber's classification of finite size systems |9j]). From this identification, we can derive many characteristics 
of the phase transition. First, consider a thermodynamic quantity for an infinite system in 3D, which varies as 
Poo(T) ~ Caot^P , where t = (T — Tc)/Tc with Tc the transition temperature for an infinite system and p an appropriate 
critical exponent. For a slab of thickness Lb, a general finite size scaling ansatz dictates that 

PlAT) -- L-BQ{L%i) (20) 

as Lb — > 00, t — > with 6 = l/v. The exponent uj is determined by requiring bulk behavior in the limit Lb -^ 00; 
this condition gives lo = pjv. The transition temperature for a finite B is shifted 

(T, - T,{Lb))/T, ^ Lg^ (21) 

and the shift exponent A is generally equal to l/v, as has been discussed for the superfluid transition in bulk "^He of 
finite thickness by Ambegaokar et al p5[ . For our purposes, it is sufficiently accurate to take v ~ 2/3, which then 
agrees very well with our numerical result (Eq. 03). Eq. Efl needs to be modified for a quantity with a logarithmic 
divergence; it becomes Poa{T) ~ Coo Ini as i -^ 0, one modifies the ansatz Q so that we have Plb {T) — Pl^ (Tq) ~ 
Q{L^gi) — Q{L^gto), where Tq is some non-critical temperature. For i ^ at a fixed Lb, we obtain for such a variable 

Pls{Tc{Lb)) ^ -CooOlnLB (22) 
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where we have assumed that Q{z) = C(l) for z -^ 0. This prediction is in good agreement with the calculated 
maximum height of specific heat peak(Eq. |l2|), which for the 3D XY model, has a weak divergence with a ~ 0. 
Similar results have been discussed for the superfluid transition in He II of finite thickness by Ambegaokar et al [^ . 
As the field increases, one may eventually reach the limit Lb/(,xy{T) < 1, at which the transition at Tc{Lb) will 
crossover to a 2D KT universality class and we expect the merging of the two transitions, Ti = Tm- In our model, we 
believe this happens for a value of / between 1/6 and 1/18. 

IX. OTHER RECENT SIMULATIONS 

Towards the completion of this work, we became aware of some of the more recent studies based on similar models. 
We briefly discuss them in comparison with our main results and interpretations. To avoid confusion, we use our own 
conventions for the flux density given in terms of the frustration / defined earlier and introduce the anisotropy factor 
r^ = Jxy/Jz- For an isotropic system, F = 1 while for YBa2Cu307_5 it is ^ 0(10) and for Bi2Sr2CaCu208, it is 
>> 0(100). We also use the same notation Tg for the temperature where jzz drops to zero. Some researchers opted 
to use Tz. 



Nguyen and Sudb0 1 35 1 have extended their earlier work on the anisotropic London loop model [ p4| . Their numerical 
results in both the vortex structure factor and jzz for F = 1 with / = 1/32 follow a pattern qualitatively similar 
to our main results for F = 1, / = 1/24 [See Figure 6 of |33]] as well as those of Li and Teitcl |23,Effl. By looking 
at the dependence of Ti{Nz) on the thickness of the system 16 < A^^ < 96, and linearly extrapolating the finite size 
effect, they conclude that T^ = Tm in the thermodynamic limit (A^^ — > oo). Their argument is based on the following 
observations: 1) T^ decreases with an approximately linear dependence on increasing A^^, Ti{Nz+5Nz) ^ Te{Nz)—c-dNz 
with a positive number c. 2) In the thermodynamic limit, below the melting transition (T ^ Tm), the energy scale for 
the interlayer phase fluctuation T* ~ ^'^Jz diverges due to long range lattice order. Therefore, the linear progression 
can not continue below Tm- From these, they conclude that Tg — > Tm in the thermodynamic limit. 

We agree with the validity of the second assumption on general grounds. However, this does not exclude the 
possible existence of an intermediate phase in which the interlayer fluctuation may be suppressed due to the quasi- 
long range phase correlations in a line liquid "phase" . With this possibility open, T„j is only a lower bound for the 
Ti. It should also be noted that "f^z does not show a significant dependence on size in the region where 0.8 < ^zz < 1, 
in temperatures Tm < T < l.STm- It is only at higher temperatures, near where jzz —^ 0, that the linear dependence 
of the shift in T^ on Nz is observable. In other words, the size dependence of ^zz is not trivial as temperature 
varies and one should not expect the same size dependence be uniformly applied over the whole temperature range 

Tm<T< TxY. 

Furthermore, the linear dependence of shift in Tg on Nz in the region where "fzz ^ is anticipated on more general 
grounds. In the London limit, "fzz measured in the simulation under periodic boundary conditions is 
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lzz{qx = n/^NxNy) ^ —- I - < Uyiq^ = 7r/Nx)ny{q + x = -Tr/N^) > (23) 
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where ny{qx — ti /Nx) is the Fourier component of the vorticity vector field lying along y direction. Position of T^ 
is governed by the condition that the vortex fluctuations make the factor in the bracket vanish. Let us assume that 
there is a characteristic number of layers A^* for which the true thermodynamic transition is realized at T^ . For a size 
Nz smaller than A^* , A^^ = A'^ + ^A'^ {5Nz < 0) , linearization of the condition gives 

Ti{Nz)^T*,+g{Tl)5Nz (24) 

with g{T) — F^9[< ny{'K /Nx)ny{—Tr / N^) > /T]/dT. What is the temperature dependence of < UyUy >? Near r„j, 
where thermally activated vortex loops begin to appear, it is dominated by the vortex loop fugacity factor and is 
steeply increasing function of T following an 5-shaped curve. However, near the foot of ^zz, where our Tg is located, 
it is numerically observed that it has reached the plateau and the temperature dependence is dominated by the 1/T 
factor. It is also consistent with the interpretation of Te in terms of XY-type unbinding transition. Since g(Ti) < 
and SNz < in this region, we then reach the conclusion that Ti{Nz) linearly increases away from T^ with c\N* — Nz\ 
as Nz decreases from the asymptotic (thermodynamic) limit. Note that we do not assume T^ = Tm in reaching this 
conclusion. If one should follow Nguyen and Sudb0 and extrapolate the observation of the linear dependence in 
the limited range of the system size and conclude that T^{B) = Tm{B), it also follows that we have a paradoxical 
consequence of predicting Txy —^ for the zero field, based on the similar bahavior numerically observed for jaa for 
/ = 0. ^ ^ 

Following the convention of Koshelev, we employ the scaled field variable h — T^ f which characterizes the ther- 
modynamics of the mixed state as long as A ^ oo. As we pointed out earlier |Q, / = 1/6 with F^ = 1 (ft, = 1/6) 
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represents a situation where the interlayer decouphng sets in right at the melting transition. From our numerical 
results with h varying from 1/6 to 0, we believe that there is a universal crossover value of scaled density h between 
1/6 and 1/18 which separates a low field regime from the high field regime for which T™ — Te ~ Tc{B). Koshelev |97| ] 
made a numerical observation that Ti — > T„j. However, it is again made for F^ = Jxvl Jz — 36, / = 1/36, i. e. h = \, 
equivalent to an extremely dense limit. Recent calculations by Hu and Tachiki pm based on / = 1/25, F^ — 10 in 
a larger system size of 50 x 50 x 40 falls into the same category with h = 10/25 ^ 1/6. For this dense limit, they 
observe that jzz drops sharply to zero at the melting transition, as wc had observed earlier for the h ~ 1/6 case in 
the stacked triangular XY model [ p3[ . 

The technical difficulty of simulating the very low field limit {h < 1/18) with large number of field induced vortex 
lines remains largely unsurmounted. To minimize the artificial grid pinning effect, one is required to choose a fairly 
small value of /( < 1/32 for the square grid, < 1/16 for the triangular grid ). However, choice of a large anisotropy 
factor to mimic HTSC then tends to push the model into the high field regime as pointed out earlier. To access 
the truly low field regime(i3 <C 1 tesla for YBa2Cu307_5 , ^ 200 gauss for Bi2SrCa2Cu208 ), we had to employ 
an isotropic model for a practical value of / < 1/18. Summarizing this section, we make educated guess on the 
thermodynamics of extremely low field limit, where the field- induced vortex degrees of freedom are no longer viable 
p9[ . Here, the transition at Te takes over the first order melting transition as the 5- A^ transition which possibly 
belongs to the universality class of zero field transition of XY film with thickness as- The phase for Tm < T < Ti 
may still be considered superconducting in the sense that a superconducting path, however narrow, may exist across 
a macroscopic distance along B, which defines a tiny, but finite critical current. This may be enhanced by collective 
pinning of the single lines, but resistance measured with a large driving will show a non-linear IV characteristics and 
hysteresis as current has to distribute itself among the superconducting paths and normal channels. 

The field induced lines in the low field limit, if we take the world- line analogy, are equivalent to extremely massive 
bosons and eventually drop out of the thermodynamics. They become localized charges which tend to polarize the 
underlying vacuum and induce a dielectric breakdown as the XY medium become more and more polarizable as T 
increases toward Tg . It is somewhat similar to a metal-insulator transition in a narrow band-gap semiconductor in 
which local field gradient (due to the field induced vortices) and shrinking band gap (as T -^ Ti) conspire to a massive 
generation of screening dipolc pairs (vortex loops) leading to a metallic state. 

X. CONCLUSIONS 

To summarize, we confirmed the existence of a single first order melting transition at high vortex densities. We also 
showed that upon melting the local vortex density increases due to the screening effect of thermally generated vortex 
loops. More significantly, however, we observed that, in the absence of disorder, destruction of phase coherence in a 
superconductor may proceed by two separate transitions at low magnetic fields: a quasi-long range phase coherence 
parallel to the field disappears at a temperature T^ higher than T„i at which the lattice periodicity disappears and 
true long range phase coherence is lost. In this low-field regime, the lattice first melts into a liquid of lines with a 
finite entanglement length along the applied field. These lines eventually disappear through increasing entanglement, 
and through their interaction with thermally induced vortex and antivortex loops. While the melting transition is 
best characterized by the disappearance of Bragg peaks for the vortex lines and a delta function peak in the specific 
heat, there is a narrow region above Tm where we observe dramatic changes in dissipation tensor which coincide with 
jump in the local vortex density and disappearance of the longitudinal phase rigidity, 7^^ — 0. Instead of being a 
gradual crossover, we propose that a possible transition at T^ ^ Tm sets in at low densities. It tracks the broad peak 
in specific heat as B increases, obeying the behavior of 3D zero field XY system confined to a semi-infinite slab of 
finite thickness Lb ^ B^^l^ . It can alternatively described in terms of appearence of connected vortex tangle which 
effectively leads to decoupling of neighbouring layers. Within this picture, origins of several puzzling and conflicting 
anomalies recently obtained on Bi2SrCa2Cu208and YBa2Cu307_5may be understood. 
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(a) 




Lz 



Lz/2 




-Vab 



FIG. 1. Geometry for the dynamic calculations described in the text, (a) To probe the shear modulus, current is injected 
uniformly into plane Ei„ and extracted uniformly from T,out- The arrows indicate directions of the Lorentz forces acting on 
the lines in the two half volumes (in opposite directions, because of the periodic boundary conditions) . (b) To probe the c-axis 
resistance, the currents are injected and extracted uniformly from the two planes indicated; the voltage drop between the planes 
is measured as in (a). 
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FIG. 2. Two illustrative phase configurations, one with net vorticity piercing the sample parallel to the xy-plane (left panel), 
and one with zero net vorticity but containing a bound vortex loop (right panel). 
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FIG. 3. 7^^ from Monte Carlo with 10,000 {T/J < 2.7) to 50,000 (T/ > 2.7) MC sweeps. The hne represents a cal- 
culation based on the harmonic self-consistent approximation. The inset shows the finite size scaling analysis to locate 
TxY = 3.04 ± 0.02 J. 
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FIG. 4. Size distribution of connected vortex segments for / = 0. The insets show typical vortex configurations for T/J = 2.7 
and 3.1. 
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FIG. 5. Dissipation (dV/dl = [V{I = 0.083) - V{I = 0.043)]/0.04) across the XY transition, Txy = 3.04J. A uniform 
current of I /Ic per grain is flown through yz planes. The inset shows p = V/I at I — 0.083/c. 
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FIG. 6. Average number of vortex segments in equilibrium(Monte Carlo) and with various bias current densities(RSJJ 
dynamics). The equihbrium density very closely follows an activated form with U = 16 J for T < Txy- 
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FIG. 7. Normalized Bragg intensity(open circles) and 72^ (filled circles) for / — 1/6. The solid line is guide for the eyes. The 
insets show the real space density correlation < nz{r,z)nz{0,0) > for the local z-vorticity taken over 40,000 MC sweeps. 
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FIG. 8. Typical flux line configurations at for two temperatures {T/J = 1.1 and 1.3) spanning the melting temperature at 
/ — 1/6, plotted for an 18 x 18 x 18 grid. The upper panels are top views. 
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FIG. 9. Specific heat per grain for various frustrations. For / = 1/6, Cv at T = 1.2 shows a clear divergent behavior, apart 
from other dilute cases. 
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FIG. 10. Dependence of peak height of Cv on the magnetic length Lb- The hne is guide for the eyes. 
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FIG. 11. Shift of peak in Cv from the / = position(T^j,) vs. the magnetic length Lb- The line(~ Lg ) is guide for the 
eyes. 
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crosses: 5 x 10'' MC sweeps) for / — 1/24. The solid lines are guides for the eyes. The inset shows the dependence of < jzz > 
on the accumulation time, as discussed in the text. 
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FIG. 13. Phase correlation function c(0; z) for various system sizes and boundary conditions. 
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FIG. 14. Longitudinal phase correlation length ^^o determined from fitting c(0; z) to the form aexp[— z/ari^o] for z < £,x{T) 
as discussed in the text. 
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FIG. 15. Real space top-to-bottom density correlation function < nz{r,z = 12)nz{0,0) > for T/J — 1.5, 1.6, 1.8 and 2.0 for 
/ = 1/24 in a 24 X 24 X 24 grid. Melting occurs around Tm/J ~ 1.55 while line-like correlations (white spot at the center) 
vanish around T/J = 2.0 ~ Te/J. 50,000 MC steps for each temperature. 
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FIG. 16. Vortex configurations for temperatures T in tlie range Tm < T < Ti. Tlie left column shows the portion of vortex 
fluctuations which form bound loops. The center column shows mainly those field-induced vortex lines which are not entangled, 
while the right column shows the largest cluster of entangled lines (as defined in the text). 
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FIG. 17. The distribution of transverse vortex lengths l^y projected onto the xy-plane l^y for three sets of temperatures in 
the (a) lattice, (b) line liquid, and (c) tangled vortex web states. 
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FIG. 18. Maximum projected lateral vortex length Ixy, given in terms ol the lateral box dimension (=48), and normalized 
per layer in a 24 x 24 x 24 grid with / = 1/24. 
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FIG. 19. Effect of switcliing connections among entangled flux lines in a torus geometry [panel (a)] and in an infinite plane 
with open boundary conditions. For the case of only two lines [(b)], note that it is necessary to make a long excursion spanning 
the whole plane to change the global winding number. The latter may be easily induced in a dense environment through 
collective occurrence of local reconnections [panel (c)]. 
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FIG. 20. Calculated bulk in-plane resistance vs. T for / = 1/24. Currents of 0.083 — 2.83/c per grain were injected uniformly 
into an yz plane. 
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FIG. 21. Calculated c— axis resistance (arbitrary units) vs. T for / = 1/24. A current of 0. 0831c per grain is injected 
uniformly into a plane. The inset shows the same results, but with temperature rescaled to model a hypothetical isotropic 
high-Tc superconductor as in the previous Figure. 
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FIG. 22. Example of the local Voronoi cell occupied by a vortex. 
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FIG. 23. Example of the vortex configuration in a xy-plane at 2 = 6 in a 48-48-12 system with / — 1/24 at T(. The black 
dots are field induced vortices, gray dots connected by lines are bound dipoles identified. The right panel shows the Delaunay 
triangulation applied to the field induced vortices only. 
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FIG. 24. Entropy S (upper panel) and normalized vortex density (lower panel) for / — 1/24 and 1/6 in 24 x 24 x 12 and 
48 X 48 X 12 systems with open boundary conditions. The error bars denote the rms deviations from layer to layer. 
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FIG. 25. Dependence of "broad" peak maximum in specific heat and their shift in temperature on the magnetic length Lb 
taken from SchiUing et al M. 
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